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ABSTRACT 


A generalized  critical  point  is 
characterized  by  totally  non-linear 
dynamics.  ^-Wo.  formulate  the  determin- 
istic and  stochastic  theory  of  relaxation 
at  such  a point.  Canonical  problems  are 
used  to  motivate  the  general  solutions. 

In  the  deterministic  theory,  we  ohow-  tbat- 
at  the  critical  point  certain  modes  have 
polynomial  (rather  than  exponential) 
growth  or  decay.  The  stochastic  relaxa- 
tion rates  can  be  calculated  in  terms  of 
various  incomplete  special  functions. 

Three  examples  are  considered.  First, 
a substrate  inhibited  reaction  (marginal 
type  dynamical  system) . Second,  the 
relaxation  of  a mea^n  ^ield  ferromagnet. 

We  obtain  g.  result,  that  generalizes  the 
work  of  Griffiths  et  al.  Third,  to-got* ^ 
sider  the  relaxation  of  a critical 
harmonic  oscillator.  , 
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SECTION  X 
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INTRODUCTION:  "CRITICAL  SLOWING  DOWN" 


Thermodynamic  and  kinetic  generalized  critical  points  are 
characterized  by  totally  non-linear  dynamics.  Such  non-linear 
dynamics  lead  to  many  interesting  phenomena,  e.g.,  "anomalous" 
fluctuations  (treated  in  (1))  and  the  "slowing"  down  of  the  decay 
of  a perturbation.  To  illustrate  the  latter  effect,  consider 
the  kinetic  equation: 

x = b(x,  a)  xCR1  cCR^  (1.1) 

for  which  the  origin  is  assumed  to  be  a steady  state:  b(0,  a)  = 0. 
Suppose  that  the  system  is  perturbed  to  a value  x = x^.  A well 
defined  problem  is  to  calculate  the  time  that  the  system  takes  to 
reach  6(0  < 6<<  xQ)  from  Xq . If  xQ  is  "small,"  then  a 
natural  approach  involves  approximating  (1.1)  by 

* 2 

x = b'(0,  a)x  + 0(x  ) x(0)  = x^.  (1.2) 

Now  we  assume  that  b'(0,  a)  < 0 , so  thatthe  perturbation  decays. 

The  time  that  the  system  takes  to  reach  x = 6 is  easily  calculated 
to  be 


i 


mm 


However,  suppose  that  for  some  critical  value  of  a = a^, 
b'(0,  ac)  = 0.  Then  (0,  ac)  is  a "generalized  critical"  point: 
the  dynamics  at  a = ac  are  totally  non-linear.  Equation  (1.3) 
yields  the  physically  ridiculous  result  t^  = 00  . Furthermore, 

(1.2)  becomes  x = 0.  Both  of  these  difficulties  are  due  to  im- 
proper linearization  procedures,  and  not  any  physical  divergences. 

In  fact,  the  decay  of  the  perturbation  is  algebraic  in  time,  with 
the  exact  form  determined  by  the  nature  of  the  singularity  at  (0,  ac) . 
Such  simple  problems  and  the  canonical  bifurcations  are  considered 
in  section  2.  The  points  essential  to  the  understanding  of  critical 
relaxation  phenomena  can  be  gained  by  study  of  one-dimensional 
systems . 

If  fluctuations  are  not  included,  a steady  state  can  not  be 
attained  in  finite  time.  Since  the  deterministic  forces  vanish  as 
a steady  state  or  equilibrium  is  approached,  the  ratio  of  fluctuation 
intensity  to  deterministic  dynamics  grows.  Thus,  the  proper  theory 
of  relaxation  must  be  a stochastic  one.  The  deterministic  kinetic 
equation  is  modified  by  a Langevin  approach.  We  use  the  diffusion 
approximation  to  treat  the  stochastic  kinetic  equation.  In  particu- 
lar, we  derive  a diffusion  equation  for  T(x'|x),  the  expected  time 
to  reach  x',  starting  at  x and  conditioned  on  the  fact  that  the 
process  reaches  x'.  We  analyze  the  one-dimensional  equations  fully 
and  obtain  certain  special  functions,  which  are  genearlized  in 
section  4 to  the  solution  of  multi-dimensional  problems.  In  sections 
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5-7,  we  consider  three  applications  of  the  theory.  In  section  5, 
relaxation  from  a steady  state  of  marginal  stability  in  a substrate 
inhibited  reaction  is  considered.  In  section  6,  we  consider  relaxa- 
tion of  a mean  field  ferromagnet.  Our  results  complement  and  extend 
the  results  of  Griffiths  et  al  (2).  Finally,  in  section  7,  we 
discuss  relaxation  phenomena  in  the  critical  harmonic  oscillator 


The  origin  is  assumed  to  be  a steady  state  of  (2.1). 

A.  Normal  Type 

The  steady  state  is  of  the  normal  type  if  b'  (0,  a)  f 0.  It 
is  stable  if  b'(0,  a)  < 0 and  unstable  if  b'(0,  a)  > 0.  In  the 
vicinity  of  the  origin,  (2.1)  can  be  replaced  by 

x = b' (0,  a)  x x(0)  = xQ  (2.2) 


As  mentioned  in  the  introduction,  the 
to  reach  x = 6,  starting  at  x = Xq 

vr 
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|b  (0,  a) 
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time  that  the  system  takes 
is 


(2.3) 
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B.  Marginal  Type 

The  steady  state  is  of  the  marginal  type  if  a(ER^  and  for 

a value  a = a we  have 
c 


b (0,  ac)  = b' (0,  ac)  = 0 , 


b (0,  ac)  t 0 


There  exists  a local  change  of  coordinates  (e.g.,  (5),  (6),  or 

Appendix  A here)  so  that  for  a near  ac  , x near  the  origin 
equation  (2.1)  becomes 


y = y - 8 (a) 


y(0)  = y0(xQ) 


(2.4) 


(2.5) 


In  (2.5),  8(a)  is  a regular  function  of  a and  8(ac)  = 0.  We 

call  a = ac  the  marginal  bifurcation  point.  The  flow  of  (2.5)  is 

sketched  in  figure  1.  The  bifurcation  picture  is  shown  in  rigure 

2.  The  marginal  case  was  considered  briefly  by  Kubo  et  al  (4) 
and  Nitzan  et  al  (7) . 

Now  suppose  that  8 > 0 and  -/fT  < yQ  < SfT  . One  can 

calculate  the  time  that  it  takes  to  reach  y^  =-/8  +6.  We  obtain 


1 2/F  | \y1  + /8)  (yQ  + /F) 


(2.6) 


For  small 


we  have 


y - /T  1 - /F/y,  2 

— = = (1  - /8/y,)  + 0(8) 

yj  + /8~  1 + /F/yj 


(2.7) 
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FIG.  1: 


DYNAMICS  OF  THE  CANONICAL  SYSTEMS 

(a)  MARGINAL  TYPE 

(b)  CRITICAL  TYPE 

(c)  HOPF  TYPE 
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FIG.  2:  (Continued) 


Expanding  the  logarithms  in  (2.6)  gives 


t ~ — + 0(B) 

*1  y0  yl 


(2.8) 


so  that  t remains  finite  as  B-0 • Clearly,  this  result  would 
Y1 

not  be  obtained  had  we  used  the  linearized  version  of  (2.5): 


- - 2 /F  ( y + /B") 


(2.9) 


In  another  possible  situation  B=0. Suppose  that  yQ  < 0 . 
The  time  to  reach  6 < 0 from  yQ  (yQ  < 6)  is  (exactly) 


1 1 

t = -T  + 

6 yn 


The  point  of  inportance  is  that  (2.8,  2.9)  yield  algebraic  forms 
for  the  relaxation  time,  whereas  (2.3)  yields  a logarithmic  time 
(i.e.,  algebraic  versus  exponential  relaxation). 


C.  Critical  Type 

2 

A steady  state  is  of  the  critical  type  if  a€R  and  for 


(2.10) 


b (0,  « ) = b (0,  a ) = b (0,  a ) = 0 

c c c 


b' ' ' (0,  a ) t 0 
c 


(2.11) 


The  canonical  form  of  the  dynamics,  for  a near  ac  and  y 

I I I 

near  0 is  (for  b < 0) 


y = -y3  + s1(a)y  + B2(a) 

y(0)  = yQ(x0) 


(2.12) 


In  (2.12),  g(a)  is  a regular  function  of  a and  6(ac)  = 0. 

We  call  a = ac  the  critical  bifurcation  point.  The  flow  of  (2.12)  is 
sketched  in  figure  1.  The  bifurcation  picture  is  shown  in  figure 
2 . When  a = ac , we  have 


y = -y' 


(2.13) 


so  that  the  origin  is  very  weakly  attracting.  The  time  that  the 
system  takes  to  reach  y = 6 from  y = > 6 is 


1_  f i _ 1_" 

2 T .2 

l*i  6 J 


As  in  the  marginal  case,  we  obtain  an  algebraic,  rather  than 
exponential,  decay  rate. 


(2.14) 


A steady  state  is  of  the  Hopf  type  if  a€R  and  when  a = a 


b(0,  a)  = b (0,  a)  = b (0,  a)  = 0 
c c c 

b'  ' ' (0,  ac)  * o 


(2.15) 
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The  canonical  dynamics  (8)  in  this  case  (for  b'1' 


< 0)  are 


where  6 = 8(a)  is  a regular  function  of  a and  6 (a  ) = 0.  The 


flow  of  (2.16)  is  sketched  in  figure  1.  The  bifurcation  picture  is 
sketched  in  figure  2.  It  is  important  to  note  the  difference  be- 


tween Hopf  and  critical  cases  (i.e.  the  number  of  parameters) 


MULTI-DIMENSIONAL  THEORY:  CANONICAL  FORMS 


We  now  consider 


we  assume 


assume  that  there  are  n distinct  eigenvalues  and  eigenvectors. 
Let  k+,  k_ , k^  denote  the  number  of  eigenvalues  with  real  part 
positive,  negative,  and  zero,  respectively.  The  dynamical  systems 


are  classified  as  follows 


A.  Normal  Case 


It  is  well  known  that  (2.17)  can  be  replaced 


by  a change  of  variables  x-*-y  so  that 


0,  then  the  origin  is  stable.  If  k > 0,  then  R can 


be  divided  into  two  sub-spaces:  an  expanding  part  (Wg) , and  a 
contracting  part  ( ) (figure  3),  with  dim  W^  + dim  W^  = n. 


B.  Marginal  Case 

We  now  let  aGR^  vary.  Then  the  eigenvalues  of  B are 
functions  of  a:  A^  = A^(a).  When  a = ac  we  assume  that, 

1)  All  eigenvalues  are  real.  Exactly  one  eigenvalue 
Ao(ac)  = 0.  There  are  k negative  eigenvalues,  A^,  ...,  X^  and 


n - 1 - k positive  eigenvalues  A,  . , . . . , X 


n-1 . 


2)  There  are  enough  eigenvectors.  Let  Z denote  the 
eigenvector  corresponding  to  An.  Then  from  (2.17),  we  obtain, 


Z = b (y  (Z)  , a) 


The  marginal  type  steady  state  is  characterized  by 


b (0 , ac)  = bz (0 , ac)  = 0 


bZZ(0'ac)  * 0 


In  appendix  A,  we  show  that  (2.17)  can  be  put  into  the  form 


(2.19) 


(2.20) 


y1  = xiy1  + 0 (y2 ) 


y°  = (y0)2  ± 6(a) 


i^nn-i 

y €R 
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FIG.  3:  MULTI-DIMENSIONAL  PHASE  SPACES. 

(a)  NORMAL  TYPE 

(b)  MARGINAL  TYPE:  DOUBLE  ARROWS  INDICATE  EXPONENTIAL 
GROWTH/DECAY;  SINGLE  ARROWS  INDICATE  POLYNOMIAL  GROWTH/DECAY 
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Our  result  is  approximate,  whereas  Arnold  (5)  and  Shashai tshvi li  (6) 
show  there  exist  transformations  which  eliminate  the  higher  terms. 

The  construction  in  appendix  A is  useful,  however,  in  that  it  gives 
explicit  ways  of  calculating  the  y and  B(a)  . When  a = ac  , 

B(ac)  = 0.  The  phase  space  Rn  is  now  decomposed  into  a direct 
product 

Rn  = W0  + We  + Wc  (2.22) 

where  Wn  is  the  manifold  corresponding  to  X„  and  W , W 
0 r 0 e c 

are  the  expanding  and  contracting  sub-spaces  respectively. 

The  assumption  that  all  eigenvalues  of  B were  real  affected  the 
form  of  the  canonical  equations.  Complex  eigenvalues  are  explicitly 
treated  in  the  Hopf  case. 

C.  Critical  Type 

2 

In  this  case,  a€R  . The  eigenvalues  of  B are  still 
denoted  by  X (a) . We  assume: 

1)  When  a = ac  there  is  one  zero  eigenvalue,  Xg, 
k negative  eingenvalues  and  n - k - 1 positive  eigenvalues. 

All  eigenvalues  are  real. 

2)  There  are  enough  eigenvectors.  Let  Z be  the 
eigenvector  belonging  to  Xg(a)  ' Then  from  (2.17),  we  obtain 

Z = b(y(Z)  , a)  (2.23) 
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The  critical  type  steady  state  is  characterized  by 


b ( 0 ,,  ac)  = bz(0,  ac ) = bzz(0,  aQ)  = 0 


bzzzf0'  “c1  * 0 


In  appendix  A,  we  show  that  the  canonical  dynamics  are 


(2.24) 


y1  = X-y1  + 0 (y2) 


(2.25) 


Z = ±(Z3)  - B1(a)Z  - B2(a)  + 0 ( Z4 ) 


In  (2.25),  B(a)  is  a regular  function  that  vanishes  when  a = a 
The  ± sign  in  (2.25)  corresponds  to  the  sign  of  b„__(0,  a ). 

Arnold  and  Shoshai tshvi li  state  a theorem  in  which  the  higher  order 
terms  are  eliminated.  As  remarked  above,  the  constructions  in 
appendix  A are  useful  for  applications.  The  decomposition  of  the 
phase  space  Rn  is  sketched  in  figure  4. 


D.  Hopf  Type 

In  the  Hopf  case,  aeR^"  and  some  of  the  eigenvalues  are  complex. 
When  ex  = one  eigenvalue,  Aq(cx)  is  a pure  imaginary  with 


35-  ReAo«» 


(2.26) 


I 


IVC 


t,<Oj  f4<0 


DECOMPOSITION  OF  PHASE  SPACE  IN  THE  CRITICAL  CASE. 

DOUBLE  ARROWS  INDICATE  EXPONENTIAL  GROWTH/DECAY.  SINGLE 
ARROWS  INDICATE  POLYNOMIAL  GROWTH/DECAY. 


r 


Thus,  as  a crosses  a^,  a pair  of  eigenvalues  closses  from  the 

left  half  plane  into  the  right  half  plane. 

i fi 

Let  x = re  . Fenichel  (1975)  (also  see  Arnold  (1972))  has 
shown  that  the  canonical  dynamics  are 

r = ±(bxr3  - hYjT)  (2.27) 

0 = X2  + b2r2  + rfY1r  (2.28) 

where  y,  is 

yi  = h Re  Va) 

A2  > 0 and  b^,b2  ^ 0.  The  function  n = n(ot)  is  regular  and 

n (0 ) = o. 
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E.  Relaxation  Rates 

Given  an  initial  displacement  from  the  origin 


Y (0)  = {yQ,  yn_x}  (2.25) 

it  is  clear  that  the  appropriate  relaxation  (or  growth)  rate  of 
the  component  (or  mode)  can  be  explicitly  calculated  by  using 

the  canonical  forms.  The  calculations  reveal  exponential  growth 
in  Wg,  decay  in  Wc  , and  polynomial  growth  or  decay  in  WQ  (at 
bifurcation  points).  Thus,  we  have  a complete,  albeit  local, 
deterministic  theory  for  relaxation  phenomena  in  the  vicinity  of 
critical  points. 


SECTION  3 


STOCHASTIC  THEORY  OF  RELAXATION: 

FORMULATION  AND  ONE-DIMENSIONAL  RESULTS 

The  deterministic  theory  of  section  2 is  approximate  in 
that  it  ignores  fluctuations-  Since  the  deterministic  dynamics 
vanish  (to  perhaps  second  order)  at  a steady  state,  the  proper 
theory  of  relaxation  phenomena  is  a stochastic  one.  Our  theory  is 
still  phenomenological,  but  it  may  be  possible  to  connect  it  to 
underlying  statistical  physics. 

3.1  STOCHASTIC  KINETIC  EQUATION  AND  DIFFUSION  APPROXIMATION 

- 

We  replace  the  deterministic  kinetic  equation  (2.17)  by 
the  Langevin  equation: 


bX(xT) 


+ H aj(x) 


y-5  (t/r ) 


(3.1) 


In  (3.1),  x is  a small  parameter  relating  the  time  scales  of  the 

fluctuations  and  the  deterministic  dynamics,  e is  a small  parameter 

characterizing  the  intensity  of  the  fluctuations.  The  process 
✓Vi  . , _ 

yJ  is  a zero  mean,  mixing  process  (for  more  exact  assumptions, 
see  (9)).  The  field  ct^(x)  is  assumed  to  be  known,  or  qiven  by 
some  prescription. 

The  process  y(S)  in  (3.1)  has  correlations.  Hence,  our 
model  is  more  reasonable  than  "white  noise"  models.  We  let 
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yk(S)y1(0)  ds 


(3.2) 


f+O,  xT^x  , a diffusion  process.  If  Uq (x)  is  a bounded. 


measurable  function  and 


u(x,  t)  = e£uq  (x(t)  ) | x ( 0 ) = , 


(3.3) 


then  u(x,  t)  satisfies  the  backward  equation 


£a  , ,,  i , i. 

u = — «—  u.  . + (b  + ec  )u.  . 


In  (3.4),  subscripts  indicate  partial  differential  and  repeated 
indices  are  summed  from  1 to  n.  The  coefficients  a1-',  c1  are 


= o^(x)oi(x)  £ykl 


i . kl  j , . 8 i , . 

c = y <ij£(x)  — ■£  o1(x) 

3x 


(3.4) 


(3.5) 


(3.6)  N ' 


In  practice  a(x)  and  c(x)  cannot  be  calculated  from  first  prin- 
ciples, but  some  prescription  must  be  given  for  their  calculation 
(e.g. , (10) ) . 

Let  N be  a neighborhood  of  a stable  steady  state  or,  more 
generally,  a domain  in  Rn  . We  set 


u(x,  t)  = 
u(x,  t)  -*■ 


u(x,  0) 


xCN 

as  distance  from  x to  N 


x CN 


(3.7) 
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Then,  u(x,  t)  is  the  probability  that  x(t)  has  entered  N by 

/v 

time  t , given  that  x(0)  = x. 

For  stochastic  relaxation  theory,  we  are  interested  in  the 
expected  time  to  enter  N , given  x(0)  = x and  that  the  process 
enters  N: 


T (x) 


t ) dt . 


(3.8) 


Then  T(x)  satisfies 


ea 


T^j  + (b1  + ec1)T^  = -u(x) 


(3.9) 


where 


u(x) 


lim  u(x,  t) 
t-»°° 


(3.10) 


Namely,  u(x)  is  the  probability  that  the  process  eventually 
enters  N,  given  that  x(0)  = x.  The  boundary  conditions  appro- 
priate to  (3.9)  are 

T (x ) = 0 x€N 

and  a growth  condition  as  distance  (x,  N)-*-°° 


(3.11) 


3.2  EXACT  SOLUTION  AND  CANONICAL  FORMS 

When  x£  R1  , equation  (3.9)  is  an  ordinary  differential 
equation 
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(3.12) 


ea 


2 Txx  + (b(x)  + e(c))Tx  = -u(x) 


Let  N = {x}  . Then  the  solution  of  (3.12)  is 

A 


T (x)  = 


X 

S 

s 

r x 

f a 

2 f b+ec  . 

f -/  X 

2 f b+ec  . 

/ 7exp 

[-  e J a dS] 

J u(x)  exp 

— QO 

e J a 

dx  ds 
(3.13) 


Equation  (3.13)  has  a rather  conplicated  asymptotic  analysis. 

Instead  of  doing  an  asymptotic  analysis  on  (3.13),  we  return  to 
(3.12)  and  set,  for  convenience  a= 2 , u(x)  = 1,  and  c=0.  We  will 
analyze  (3.12)  and  obtain  certain  special  functions.  These  functions 
will  be  generalized  in  section  4,  for  the  solution  of  multi-dimen- 
sional problems.  Our  analysis  is  based  on  matched  asymptotic 
expansions  (e.g. , (11)). 

Away  from  the  zeros  of  b(x),  (3.12)  becomes 

b(x)Tx  = -1  (3.14) 

This  is  the  "outer"  equation. 

Near  zeros  of  b(x),  (3.14)  breaks  down.  We  need  to 

stretch  coordinates  (3.12)  to  obtain  the  appropriate  "inner" 
equations.  We  shall  analyze  (3.12)  by  using  the  canonical 
form  of  b (x) . 
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A.  Normal  Case 

In  the  normal  case,  b(x)  = ±x,  with  (+)  indicating  that 

the  origin  is  an  unstable  steady  state,  (-)  indicating  a stable 

steady  state.  Introducing  z = x//e  , (3.12)  becomes 

T ± zT  = -1 

zz  z 


(3.15) 


B.  Marginal  Case 

. 2 ~ 
In  the  marginal  case,  the  canonical  dynamics  are  b(x)  = ± (x  )-  u, 

We  introduce  the  stretched  variables 
1/3  ~ 2/3 

z = x/e  , a = a/e  , so  that  (3.12)  becomes 


T ± (z2  - a ) T = -l/e1/3 

z z z 


(3.16) 


C.  Critical  Case 


In  the  critical  case,  the  canonical  dynamics  are  b(x)  = 

3 ~ ~ 

± x + 6-^  x+$2  ' We  intr°duce  stretched  variables'  z = x/e 4 , 


61  = 6l/eJ5  ' 62  = B2/e3/4 


and  obtain  the  inner  equation 


T „ + (±z3  + B.z  + e_)T  = - 1/e*5 

ZZ  i 2 z 


(3.17) 


D.  Hopf  Case 

3 ~ 

In  the  Hopf  case,  the  canonical  dynamics  are  b(x)  = -x  + gx  . 

v ~ l 

We  introduce  the  stretched  variables  z = x/e 4 , 6 = 6/e^  and 
obtain  the  inner  equation 


T + ( — z 3 + Bz)T  = - 1/e*3 
zz  z 


(3.18) 


i 


Equations  (3.15-3.18)  define  certain  incomplete  special 
functions.  These  special  functions  will  be  used  in  the  next 
section  to  construct  asymptotic  solutions  of  multi-dimensional 
problems . 


SECTION  4 


STOCHASTIC  THEORY:  ASYMPTOTIC  RESULTS 


When  x€R  ns2  , equation  (3.9)  will  usually  not  have  exact 
solutions .Consequent ly , approximate  techniques  are  required.  The 
methods  used  here  are  closely  related  to  those  in  (10) . The  basic 
idea  is  to  generalize  the  one-dimensional  inner  solutions;  we  call 
the  method  a generalized  ray  method.  Although  the  normal  case  does 
not  represent  a "critical"  point,  we  include  it  for  completeness. 

4.1  NORMAL  CASE 

We  suppose  that  the  origin  is  a simple  steady  state  (figure  5) 
and  that  it  is  stable.  We  seek  a solution  of  (3.9)  in  the  form 

T (x)  = g(x)F((| + htxje*5  f'(^/v'7)  + k(x)  . (4.1) 

In  equation  (41.),  F(z)  is  a special  function  satisfying 

d2F  dF  , , „ 

2 Z dz  ~ 1 (4.2) 

dz 

and  the  functions  iji  (x)  , g (x)  , h(x)  , and  k(x)  are  to  be  determined. 

In  order  to  completely  analyze  the  problem,  we  assume  that  g,  h,  k 
have  expansions 

g(x)  = 2 gn(x)en  h(x)  = I hn(x)en  k(x)  = 2 ka(x)en  (4.3) 

Consequently,  the  construction  given  here  represents  the  first  term 
in  the  asymptotic  solution  of  (3.9). 

When  derivatives  are  evaluated,  (4.2)  is  used  to  replace 
F"(t|r/,/e)  by  F ' (♦/-/£)  - 1.  Then  terms  are  collected  according 

to  powers  of  e.  We  obtain: 
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-u(x)  = 


_ -h 


e [b1  i|i  . + a1-1  (g+h\|i ) F ' 

2 3 


+ e°F(bigi) 


ID 


+ e°  (b1ki  + 4>i4' jg) 

ij 


+ e^F ' ^b1hi  + 9>lij  + aJUJhi>)i  j f 

i j i j 

+ "2  h't'ij^+  ~y“  titjh  - gc1^ 

i*i*] 


+ he 

The  leading  terms  vanish  if 
i aij 

b ’I'i  + — 'I'i'l'j#  = 0 

b1gi  = 0 


i a13 

b ki  + 2 ^ ± 'l' -j g = -u(x) 


(4.4) 


(4.5) 

(4.6) 


(4.7) 


First  consider  (4.5).  Since  b1(0)  = 0 for  all  i,  we  set  \(i(0)  = 0, 

in  order  to  keep  i|i(x)  regular.  Then  (4.5)  can  be  solved  by  the 

method  of  characteristics.  We  note  that  the  transformation 
2 

$ = *si|i  converts  (4.5)  to 

ij 


, 1.  . a 
b + — 


2"  *i*j  = 0 ' 


(4.8) 


which  is  a Hamilton-Jacobi  equation  (see  also  (12)).  Then,  we  can 
solve  the  Hamilton-Jacobi  equation  in  terms  of  characteristics: 


v1  = 
- 


;i  _ -dH 

p - “1 
ax 


* = \ al3PiPj 


(4.9) 
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Starting  at  the  origin,  the  phase  plane  is  covered  with  trajectories, 
called  rays,  along  which  (or  $)  is  known.  Thus,  ^ at  any  point 

x is  known. 

Equation  (4.3)  indicates  that  g is  constant  on  deterministic 
trajectories.  Since  all  trajectories  intersect  at  the  origin,  g 
must  have  the  same  value  on  all  trajectories.  At  the  origin,  (4.7) 
becomes 


T"  Vjg 


= -u  (0 ) = 


(4.11) 


Thus 


(4.12) 


V- 

£ l 3 


We  set  MO)  =0  as  initial  data  for  (4.7). 

If  we  set  F ( 0 ) = F * ( 0 ) =0  as  initial  conditions  in  (4.2), 
then  the  leading  term  of  the  asymptotic  solution  satifies  T(0)  = 0 


The 

0 u1/2 

bV 

l 

+ «ij 

+ 2 

+ -j-  - gc  iJk  + he  = 0 
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(4.13) 


At  the  origin  b'  (0)  = ip  ( 0 ) = 0,  so  that  (4.13)  becomes 
1 


h (0 ) = — 
a 


|-a1J  , x i, 

\ ^ — W . . g + C lb  . q 

3,*,  \ 2 


(4.14) 


J-  3 


Equation  (4.13)  can  be  solved  by  the  method  of  characteristics,  with 
initial  data  given  by  (4.14). 

Thus,  we  have  completely  constructed  the  leading  term  of  the 
asymptotic  solution  of  (3.9). 

As  a by-product  of  our  method,  we  are  able  to  approximately  solve 
the  famous  Kolmogorov  first  exit  problem,  recently  considered  by 
Matkowoky  and  Schuss  (13)  using  matched  asymptotic  expansions.  This 
problem  is  the  following:  suppose  that  the  origin  is  surrounded  by  a 
domain  D , with  boundary  3D.  Find  the  expected  time  that  the 
process  takes  to  hit  the  boundary  (i.e.  the  mean  exit  time  from  D) 
(Fig.  5b)  from  x. 
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We  follow  the  arguments  leading  to  equations  (4.1  - 4.11), 
except  that  the  initial  data  for  F,F'  and  k(x)  change.  We  set 
k(x)  h 0 on  3D.  We  distinguish  two  cases: 

i)  The  boundary  dD  is  a contour  of  i|/  (or  $)  say,  on 

dD.  Then  we  set 

F(!|id/77)  = F'(^D/Vl)  = 0 (4.15) 

when  solving  (4.2).  Then  T s 0 on  D. 

ii)  The  boundary  3D  is  not  a contour  of  f.  Let  iji^  and  ^ 

denote  the  maximum  and  minimum  values  of  | on  3D.  Then  T | 0 on 
3D,  but  it  can  be  shown  that  on  3D 

|T|s|&,(i|iI/|iII)|  (4.16) 

+ exponentially  small  terms. 

Hence,  if  \Bn{h ) | is  small,  then  |T(x)  | will  be  small  on  the 
boundary. 

4.2  MARGINAL  CASE 

In  some  senses,  the  marginal  case  has  the  least  interesting 
dynamics.  The  dynamical  problem  we  consider  here  is  sketched  on 
figure  6.  When  the  deterministic  system  has  two  nodes  (QQ,  Q-^) 

and  one  saddle  (S) , even  if  the  process  starts  near  Q , it  will 
• eventually  reach  Q^,  due  to  the  proximity  of  Qq  and  S.  The 

proper  question  in  the  stochastic  theory  involves  the  time  to  cross 
some  given  curve  R.  We  note  that  such  a time  is  infinite  in  the 
deterministic  case,  if  the  phase  point  starts  on  or  above  S. 
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• tf* 


We  seek  a solution  of  (3.9)  of  the  form 

T (x)  = g(x)B(f/e1/3,3/e2/3,l/e1/3, y2) 

+ h(x)e2/3B,($/e2/3,l/e1/3,Y2  ) + k(x) 

In  (4.17),  B(z,  “,X1,X2)  satisfies 


(4.17) 


d2B  = ,„2  . dB  . , „ 

2 (Z  a)  dz  X1  X2 
a z 


(4.18) 


and  g(x),  h(x),  k(x),  t|/(x)  and  the  parameters  a»Y2  are  to  be  deter- 
mined. We  proceed  as  in  Section  4.1.  Instead  of  equations  (4.5-7) 
we  obtain 

i aij  2 

b (*  - 30)  = 0 

b^g.  = 0 

3 L 

i aij 

b ki  - -y  ^(i.gfl  + "Y 2 ’■V ^ = 

In  (4.17),  we  have  set  g = Eg,  £ . 

K 

2 

We  set  tj/  = at  QQ  and  at  S.In  particular  f(QQ)  = + v|3q 

and  \)i(S  ) = “./go  • The  value  of  g^  can  be  determined  by  an  itera- 
tive procedure  (10) . We  pick  an  initial  value  of  gQ  = and 

solve  (4.19)  by  the  method  of  characteristics,  starting  at  Qq,  where 


(4.19) 


(4.20) 


(4.21) 


* = 


Some  rays  will  approach  S . As  a ray  approaches  S , 


ij»  should  approach  -,/gg0)  . If  it  does  not,  then  the  g^0)  must  be 
replaced  by  a second  iterate  g*1*.  The  method  of  false  position 


can  be  used  to  calculate  iterates  of  3 


This  procedure  can  be  re- 


0 

peated  until  30  is  known  to  any  desired  accuracy.  (In  (10),  we 
present  a discussion  of  this  calculation  in  more  detail.) 

Equation  (4.20)  indicates  that  g is  a constant.  At  Qq  and 
Q^,  which  we  denote  generically  by  P , we  have,  from  (4.21)  : 

g ( 1 + Y 2 ’MP)  ) = -u(P)  (4.23) 


— ip.ip. 

2 ri  j 


These  are  two  equations  for  the  unknowns  g and  y2-  We  set  k = 0 
on  R and  assume  that  R is  a level  curve  of  ip , with  tp  = <pR  on  R. 
Then  we  set 


B(<J>R/e1/3, 3, 1/e1/3,Y2)  = B1  ( ^R/e1/3 , 3 , l/e1/3  , Y2)  = 0 (4.21) 

With  these  choices,  T(x)  =0  if  x R, 

At  the  bifurcation  point  n = 0 (the  marginal  bifurcation) 

Qq  and  coalesce.  Then  3q  = 0 , and  it  can  be  shown  that 

Y2  = 0 (10).  At  the  saddle-node  Qq/Q^,  equation  (4.23)  still  pro- 
vides one  equation  for  g: 


g 


u(P) 


13 


ip . ip  . 


(4.24) 


Elsewhere,  we  have  given  proofs  that  all  the  construction  are  regu- 
lar at  the  bifurcation  point  ((10),  appendices  D,  E) . 

In  section  5,  we  consider  an  example  of  a chemical  system  ex- 


hibiting the  marginal  bifurcation. 


4.3  CRITICAL  CASE 


We  now  consider  a system  with  three  steady  states,  Pq,  P^,  and 
P2  when  ai'a2  > When  = a2  = 0 the  three  steady  states 

coalesce  into  a critical  type  steady  state.  When  a^,a2  < 0 there  is 
only  one  real  steady  state;  it  is  assumed  to  be  stable.  If 
a^,«2  > 0,  we  surround  P2  by  a domain  N and  pose  the  following 
stochastic  relaxation  problem:  what  is  the  expected  time  to  enter  N, 
given  the  initial  position.  Clearly  there  is  an  analogous  problem 
for  a neighborhood  N of  PQ.  When  there  is  only  one  steady  state 
P,  we  surround  P by  N.  We  note  that  if  N shrinks  to  P , 
then  we  have  the  expected  time  to  "reach"  P,  conditioned  on  initial 
position.  We  also  note  that  T(x)  =0  if  xCN. 

We  seek  a solution  of  (3.9)  in  the  form 

T (x)  = g(x)Q(ii»/e1//4,n/e1//2,  B/e3//4/l/e1^2/Yi/e1'/4fY2^ 

+ h(x)e3/4Q'(^/e1/4/a/£1/2,B/e3/4,l/e1/2,Y1/6l/4,ii>2)  (4.25) 

+ k (x) 


where  Q ( z , a , 8 , » Y2 > Y3)  satisfies 


a Q . , 3 dQ 

TT  = 4 (z  - az  “ B)  dz  ' Yi 


+ y0z  + y.z 


(4.26) 


The  (+)  sign  in  (4.26)  corresponds  to  the  steady  state  P being 
stable,  the  (-)  sign  to  it  being  unstable.  We  consider  the  case  in 
which  P is  stable. 

Instead  of  (4.5),  we  obtain 


1 J o 

b1^.  + \p.ip.  (t|>  - ai|)  - 8)  = 0. 


(4.27) 


When  there  are  three  steady  states,  a and  3 are  determined  by  a 
procedure  analogous  to  the  one  in  section  4.2.  Namely,  at  the  steady 
states  we  set 

iji3-aiji-p=0.  (4.28) 

The  method  of  characteristics  is  then  used  to  determine  a and  p 
by  an  iterative  procedure.  When  the  three  steady  states  coalesce 
a = p = 0.  When  there  is  one  real  and  two  imaginary  steady  states, 
then  a,  P < 0 and  can  be  determined  by  power  series.  Such  series 
are  constructed  elsewhere  (10) . 

Instead  of  (4.7),  we  obtain 

aij  2 

b ki  + ~2~  jg  ( -1+  Y2 1If  + Y3  'l'  ) = -u(x)  (4.29) 

At  the  steady  states,  we  obtain 
aij  2 

~2~  fii|i  (-1  + Y2t  + Y3>li  ) = -u(x)  (4.30) 

When  there  are  three  real  steady  states,  we  obtain  three  equations 
for  g,  , and  y^.  When  two  steady  states  coalesce,  - 0.  We 

still  have  two  equations  for  g and  y . Finally,  when  all  three 

coalesce,  y2  = Y3  = 0 and  we  are  left  with  one  equation  for  g. 

We  obtain  an  equation  for  h(x)  that  is  analogous  to  (4.13), 
and  is  treated  in  an  analogous  fashion.  The  initial  values  of  Q 
and  Q’  in  (4.26)  are  determined  so  that  T(x)  = 0 if  x€dN. 
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4. 4 HOPF  CASE 


The  Hopf  type  dynamical  system  is  treated  in  an  identical 
fashion  to  the  marginal  and  critical  type  systems.  We  seek  a solu- 
tion of  (3.9)  in  the  form 


T(x) 


= g(x)H(.J;/e1/4,e/e1/2,l/e1/2/Y2/^l/4) 

+ VV  (4'/e1/4,l/e1/2,3/E1/2,Y2/fc.i/4)h(x) 


(4.31) 


+ k(x) 


whe 


re  H (z, P,y^/Y2)  satisfies 


— ? = ± (z3  - 3z)  - Y,  + Y02 

dt  dz 


(4.32) 


The  (+)  sign  corresponds  to  a stable  limit  cycle  and  unstable  focus, 
the  (-)  sign  corresponds  to  an  unstable  limit  cycle  and  stable  focus. 
The  analysis  proceeds  exactly  as  in  section  4.2,3. 


SECTION  5 


SUBSTRATE  INHIBITED  REACTIONS: 
A MARGINAL  TYPE  STEADY  STATE 


The  following  equations  model  a substrate  inhibited  chemical 


reaction  in  an  open  reactor  (10,14): 

'1  -1  4x^  1 

x = , -y  - .069979  x + . 25901- 

1. 5+x  +13  (x  ) 


" 2 

x = .09- 


1 2 
X X 

1 2 
1+1  Ox  x^ 


1 2 
-X  X 

1 2 

1+lOx  x 


(5.1) 


(5.2) 


1 2 

where  x and  x are  dimensionless  "concentration"  variables.  The 
steady  state  (.4359,  2.065)  is  a saddle  node,  it  is  a marginal  type 
steady  state.  The  steady  state  (1.46,  .52)  is  a stable  node.  The 
phase  portrait  is  shown  in  figure  8,  along  with  a first  exit  boun- 
dary. The  theory  on  section  4.2  applies.  We  wish  to  calculate  the 


expected  time  to  hit  R,  conditioned  on  initial  position.  Using  the 
birth  and  death  approach  to  chemical  kinetics  (15)  , can  be  modeled 
as  (10): 


where 


Ux+  bL)x 


1 2 
x x 

l+10x1x2 


, , . 1.4x 

(X-i+Pl^x  1 . 

1 1 1.5+x  +13(x  ) 


U2H2): 


= .09  + 


1 2 
x x 

1 2 

1+1  Ox  x 


(k2+p2) x" 


+ . 069979x  + .25901  + 


1 2 
x x 

1 2 

1+lOx  x 


1 2 
X x 

1 2 

1+lOx  x 


(5.3) 


(5.4) 


(5.5) 


THE  BOUNDARY  R WAS 
EXIT  TIME. 


The  parameter  e characterizes  the  intensity  of  fluctuation. 


In  table  I,  we  compare  the  theory  of  section  4 with  Monte  Carlo 
experiments  for  e = .01. 


TABLE  1 


Comparison  of  the  Theory  and  Monte  Carlo  Experiments 
in  the  Marginal  Bifurcation 


Test  Point 

T(x)  Theory 

T(x)  Experiment  (//  Trials) 

(.42,  2.06) 

60.3 

36.4  (950) 

(.38,  2.36) 

104.1 

91.2  (-400) 

/— N 

o 

fsl 

o 

rg 

66.1 

62.4  (2000) 

(.3,  1.8) 

37.7 

35.0  (1550) 

(.16,  2.4) 

119.6 

103.5  (400) 

(.7,  2.2) 

36.1 

31.4  (1750) 

(.6,  2.4) 

74.9 

b8 . 2 (800) 

SECTION  6 


KINETIC  MODEL  OF  THE  FERROMAGNET 


We  shall  give  an  analysis  of  the  mean  field  ferromagnet, 
similar  to  that  of  Griffiths  et.  al.  (2) . The  problem  is  one 
dimensional,  so  that  the  full  theory  of  section  4 is  not  needed. 
However,  this  application  illustrates  many  of  the  ideas  that  run 
through  an  analysis. 

Consider  N spins,  with  cn  = ±1,  in  a magnetic  field  H. 
Let  J be  a coupling  constant.  The  Hamiltoman  is 

“ - ts  - 1,11  E°i  - 1/2J 

l<n 


(6.1) 


We  let 


n = ^ (N  + Loi)  , 


denote  the  number  of  spins  "pointing  up."  Then  (6.1)  becomes 


(6.2) 


ft  = *(n)  = — ('2-^— — yH  ( 2n  - N) 


We  take  a mean  field  approach  and  assume  that  the  number  of  spins 


(6.3) 


pointing  up  is  really  a statistical  variable,  n(t).  The  statistical 


behavior  of  n(t)  is  described  by  transition  probabilities: 


Pr|n(x  + Sx)-n(x)  = ljn(x)  = n| 

= N ~ n exp  (<Mn  + 1)  - 4>(n))J 


(6.4) 


6x  + o (6x) 


i 


Pr 


(t  + 6x)-n(x)  = -1  n(x)  = 


n; 


n 

= - exp 


(<J>  (n  - 1)  - 4>  (n) ) 


6x  + o(6x)  , 


(6.5) 


where  8 = l/knT.  We  assume  that  the  probability  of  all  other  tran- 
sitions  is  o(<5x).  In  deriving  (6.4,5),  we  have  restated  the  argu- 
ment in  (2) . We  follow  (2)  and  introduce  a "continuous"  variable 


-v...  E°i  2n  - N 

x (t)  = 


N 


N 


(6.6) 


If  6x=  x(x  + 6x)  - x(x),  then  (6.4,5)  become 


Pr  | 6x  = 2/N  | x ( x ) = xj  = 1 -x  exp  |-8(-xJ  - - Hy)j 


Hy )>  6x 


+ o(6x) 


(6.7) 


Pr  1 6x  = -2/N  x(x)  = x|  = 1 ~ x exp  |-B(xJ  - + Hy)| 


- + Hy) } 6x 


+ o(6x) 


(6.8) 


We  set  a = J8,  f>  = 6yH  and  introduce  a macroscopic  "physical" 
time  defined  by 


‘■i 


(6.9) 


Thus,  we  construct  drift  and  diffusion  coefficients 


= (1  - x)  exp  £ax  + ^ + ~ (1  + x)  exp 


(6.10) 


-ax  + 


H 


(6.11) 


b (x) 


r 


¥** 


FIGURE  9:  MEAN  FIELD  FERROMAGNET 


37a 


and 


a (x)  = 


lim 

6t-*0 


E j(6x) 2 |x(t)  = x | 


(6.12) 


= ^ | (1  - x)  exp  (ax  + ^ + 6)  + (1  + x)  exp(-ax  + ^ - 6)| 

(6.13) 


Thus,  the  average  value  of  x(t)  evolves  according  to 

x = b(x,a,6)  = 2ea/^N  jsinh(ax  +6)  - x cosh  (ax  +6)j  , 


(6.14) 


subject  to  -1<x<1.  The  steady  states  and  true  (physical)  equilibrium 
are  solutions  of  b(x,a,6)  = 0.  Therefore,  one  obtains 


x = tanh  (ax  + 6) 


(6.15) 


Equation  (6.15)  is  usually  obtained  by  a statistical  thermodynamics 
argument  (e.g.  (16),  pg.  101). 

This  agreement  adds  support  to  our  statistical  approach.  In 
many  respects,  the  approach  used  here  is  preferable  to  the  standard 
approach.  Not  only  does  the  stochastic  approach  yield  the  equili- 
brium solution,  it  gives  dynamics  and  the  steady  states.  As  is  well 
known,  equation  (6.15)  may  have  1,  2,  or  3 solutions,  depending  upon 
the  values  of  a and  6.  In  figures  9a,  b,  we  illustrate  the 
graphical  solution  of  (6.15)  for  zero  field  (6  = 0).  When  6=0, 
Xq  and  *2  are  both  thermodynamically,  and  kinetically,  stable. 
However,  for  6^0,  one  of  xQ , X£  becomes  kinetically  stable 
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(thermodynamically  metastable)  while  the  other  is  the  true  thermo- 
dynamic (and  kinetic)  equilibrium  (Figure  9c) . The  kinetic  condi- 
tion of  criticality  is  that,  when  6=0 


b' (xx)  = b" (xx)  = 0 


(6.16) 


We  easily  obtain  a = 1 as  the  critical  value  of  a.  This  defines 
the  critical  temperature. 

Now  consider  6/0,  with  xQ  metastable  and  x2  stable. 

'Xj 

The  expected  time  to  reach  x2,  given  that  x(0)  = x satisfies 


-1  = §T 

+ bT 

(6.17) 

2 xx 

X 

T(x2)  = 0 

lim 

T (x)  <“> 

(6.18) 

X-*  -°° 

a(x)  and 

b(x)  given 

by  (6.13)  and  (6.11). 

Define  the  re- 

laxatdon  rate  from  the  metastable  to  stable  state  by 

k = - 

T(xq)  * 


(6.19) 


We  can  calculate  the  relaxation  rate  k for  all  values  of  N.  The 
method  of  Griffiths  et.al.  (2)  broke  down  for  large  N.  The  result 
given  here  will  be  valid  for  all  values  of  N.  Our  result  thus 
extends  their  analysis.  It  can  be  shown  that  the  two  results  are 
equivalent  for  small  N. 
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SECTION  7 

RELAXATION  OF  A CRITICAL  HARMONIC  OSCILLATOR 


The  application  in  section  6 did  not  use  the  theory  of  sec- 
tion 4,  but  the  one  in  this  section  does.  We  consider  a Duffing 
oscillator 


dx 


dt 


= v 


md  v 


- (-k(n)x  - a(x3)  - yv)  + J ea  ^ 


dt  ' 1 v'  - dt 
We  assume  that  k(r|C)  = 0 for  some  critical  value  of  g and  that 
k(g)^0  for  all  g.  The  mean  motion  of  the  oscillator  is  given  by 
x = v 


. -kx-gx  -yv 


n 


When  a>0  , the  origin  is  the  only  real  steady  state.  The  matrix 

(°  1 

B ■ (b  'i'lo.o  =1  v 


ZX 

m 


has  eigenvalues  and  eigenvectors 


*■  ♦ = 


. - yi^/y^-4k 
2m 


e.  = 


4k 


2m 


(7.1) 

(7.2) 


(7.3) 

(7.4) 


(7.5) 


(7.6) 


when  k £ 0 , the  origin  is  a critical  type  steady  state.  Accord- 
ing to  the  fluctuation  dissipation  theorem,  for  this  problem 


a = 2kT y p , where 
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l 


•s-  v» 


E(y(S)y(0)  )ds 


(7.7) 


Let  T(x,v)  be  the  expected  time  that  the  process  takes  to  enter  a 
small  ellipse  around  the  origin,  given  that  x(0)  = x,  ^(0)  = v. 


Then,  for  arbitrary  k , 


-1  - ^ Tvv  + VT 
m 


(kx+ax+vv) 


We  introduce  scaled  variables  by 


v. 

w m 


X = — - t 1 

Y0 


Y = Ynri(x'  ) 


t = — t* 

Y0 


VEQm 


(7.8) 


(7.9) 


'°  V E0m 


Where  Eg  is  some  reference  energy,  such  that  pkT«EQ.  Defining 
e = pkT/Eg  , we  obtain  (for  k ( rj)  5 0) 

-1  = evTv,v,  + VT'x,  - (a'  (x ' ) 3 + n,v')T’v, 

In  the  sequel,  we  drop  the  primes.  Since  the  origin  is  a critical 
type  steady  state,  the  theory  of  section  4 applies. 

The  leading  term  in  the  asymptotic  solution  of  (7.10)  is 

T (x)  ~g°  Q(f(x,v)/e1/4,,0,0,  1/e1/2,0,0) 

+k°(x)  + 0 ( e3/4 ) 


(7.10) 


(7.11) 


Equations  (4.27)  and  (4.29)  become 


3 2 3 

v\|/x  -(ax  + nv)*v  + T]t v>|f  - 0 


i u / J . x , U U V 

vkv  - (ax  + r|v)kv  - g ^ 


v 


= -1 


(7.12) 


(7.13) 


In  order  to  keep  f regular  at  (0,0),  we  set  ty  = 0 there.  In 
order  to  solve  (7.12)  by  the  method  of  characteristics,  we  need  in- 


itial data  for  \Ii  and 
rx 


If  (7.12)  is  differentiated  with  re- 


spect to  v and  evaluated  at  (0,0),  we  obtain 
*x  - = 0 at  (0 , 0) 

When  (7.12)  is  differentiated  three  times  with  respect  to  x and 
evaluated  at  (0,0),  we  obtain 
tj!*v  = a/rl 


(7.14) 


Thus  we  obtain,  at  (0,0) 
*x  “ (arl)  1/5 


* 


V 


4i/^4/s 


-i 


(7.15) 


(7.16) 


Higher  derivatives  are  evaluated  in  a similar  fashion.  Thus,  we 
can  specify  an  ellipse  around  the  origin: 


N = | (x  , v)  : i|i  ( x , v)  = 6 | 


(7.17) 


We  set  Q(6/e1/4,0,0  ,l/e1/2,0  ,0  ) = Q ’ ( <5/e1/4 , 0 , 0 , l/e1/2 , 0 , 0 ) = 0 
when  integrating  (4.26).  We  also  set  k(x,v)  =0  if  (x,v)£N. 
ht  the  origin,  (7.13)  becomes 


0 

9 = 


Tl*! 


2 3/5 

~275  n 

a 


(7.18) 
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i 


J 


which  determines  the  valup  of  gJ.  Then,  on  deterministic  trajec- 
tories we  have 


dk 


0 


0 ,2 


dt 


- = -1  + 


g Tit 


v 


with  the  initial  data  given  above.  Equation  (7.12)  can  now  be 
solved  by  the  method  of  characteristics,  so  that  the  leading  term 
in  the  asymptotic  solution  is  known. 


(7.19) 
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Appendix  : Classification  of  Two  Dimensional  Dynamical  Systems. 

The  Normal  Forms  for  Marginal  and  Critical  Type  De- 
terministic Systems 

In  this  appendix,  we  classify  planar  dynamical  systems  into 
normal,  marginal,  or  critical  types.  Our  scheme  is  a generalization 


of  the  work  of  Kubo  et.al.  (4).  We 


use  Segre ' s method  (17)  to 


derive  the  local  normal  forms  for  the  marginal  and  critical  cases. 

I.  Normal  Type  Dynamical  Systems 
Let  the  dynamical  system 

x = b(x)  x€R^  ( 

have  a steady  state  at  x = x^.  Let  \±  denote  the  eigenvalues  of 

(b^J  evaluated  at  x^.  The  system  (A.l)  is  of  the  normal  type  if 

the  real  parts  of  \±  are  non-zero.  The  steady  state  is  stable  if 
the  real  parts  of  \±  are  negative. 

According  to  the  standard  theory  of  differential  equations  (18) 
there  exists  a change  of  variables  x-*y  so  that 


(A.l) 


X+0  \/ y 


0 X-/\y 


2 2 
ayl  + b yiy2  + Cy2 


a2y2  + b2yiy2  + c2y2 


+ 0 (y  ) 

The  coefficients  a^  - c^  are  given  in  terms  of  the  second  deriva- 
tives of  b(x)  evaluated  at  xQ. 
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II.  The  Marginal  Type  Systems 
The  dynamical  system 

x = b(x,n)  x<ER2  n€R  (A. 2) 

is  assumed  to  have  the  following  behavior.  (A. 2)  is  assumed  to  have 
three  real  steady  states  for  r)>0.  We  denote  these  by  QQ(r|)#  Qj_(n) 

P • 

2 

Denote  by  the  matrix  (b^)  evaluated  at  Qq,Q,  or  p?*  (k=0,l,2). 

We  assume 

1)  For  all  r),B2  has  two  real  negative  eigenvalues 

2)  As  nl0,  the  distance  between  QQ(r])  and  Q -^  ( r))  decreases. 

When  r)  = 0 the  two  points  coalesce  and  annihilate  each  other 

3)  For  t]>0,Bq  has  two  real,  negative  eigenvalues  which  de- 
pend upon  r|  and  B j,  has  one  real  positive  and  one  real  negative 
eigenvalue.  When  r|  = 0 , B^  = B^  has  one  zero  and  real  negative 

eigenvalue.  The  eigenvector  corresponding  to  the  negative  eigen- 
value has  positive  slope. 

We  introduce  as  new  coordinates  the  eigenvectors  y^,y^  aS 
new  coordinates  so  that  (A. 2)  becomes 
i'l  = ^(Y/Ti) 

y 2 = E2(y'n) 

The  system  (A. 2)  is  of  the  marginal  type  if  the  above  condi- 
tions hold  and 


( A. 2a) 


A-2 


1 


then 


2 = Y'2  = ^0^2  + 0(t>Z  > 


zi  = y{  + syiy2  + sy i^2  + 2t^iyi 


(A. 7) 


+ 2uy2y^  + wy^n  + ry^r| 

+ 0 (nyy’ ,y2y ' ) 


Using  (A. 5)  in  (A. 7)  yields 


Z2  = X0Z2  + 0(T1'Z  > 


= c^n  + y1(c2ri  + 2tc1ri)  + Y1 

+ y^^*13  + sA.q)  + y2(rr)\  + sc^]) 


+ y2  ( 2uX  + c)  + 0 (y3  , t]y2 . T)2) 


We  choose 


(A. 8) 


2 2 3 

and  note  that  y^  = + 0 (y  ).  Thus  equation  (A. 8)  becomes 

* V2  + ° |T''Z2) 

Z[  = z2  - p ( n)  + 0(z3,nz^>T12) 


(A. 9) 


where  g(q)  = -c  ^t).  Equation  (A. 9)  is  the  local  normal  form  which 
we  desire.  It  is  a weaker  result  than  that  of  Arnol’d  (5)  or 


Shoshaitshvili  (6)  who  actually  eliminate  the  higher  order  terms. 
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III.  Critical  Type  Systems 


We  now  consider  a dynamical  system  depending  upon  two  para- 


meters 


x = b (x, q, 6) 


T] , 6€R 


(A. 10) 


We  make  the  following  assumptions: 

1)  For  some  combinations  of  r)  and  6 , equation  (A.  10)  has 
three  steady  states  PQ(q,&),  P1<ti,6)  and  P2(q,6).  When  the 
three  points  are  distinct,  we  assume  that  PQ  and  P2  are  stable 
nodes  and  that  P^  is  a saddle  point. 

2)  As  q,  6 vary,  two  of  the  points  may  coalesce  into  a point 
of  neutral  stability  (i.e.,  one  eigenvalue  of  the  linearized  equa- 
tions is  zero) . This  situation  is  equivalent  to  a marginal  type 
dynamical  system. 

3)  As  q,6-*0  from  above,  the  three  steady  states  approach  each 
other  and  coalesce  when  q = 5 = 0.  Let  B = (fc>,j)  evaluated  at 

P^.  When  q,  6>0  , we  assume  that  B has  one  real  positive  and  one 

real  negative  eigenvalue.  When  q = 6 = 0 , B has  one  real  nega- 
tive and  one  zero  eigenvalue.  When  q, 6<0  , B has  two  real  nega- 
tive eigenvalues. 

We  denote  by  y^(q,6),  y2(q»6)  the  eigenvectors  of  B. 


The  eigenvectors  are  introduced  as  new  coordinates  so 

that  (A. 10)  becomes 


i 


A- 5 


(A. 10a) 


Y1  = (y , rp  6) 

y 2 = E2  (y  / n»  6) 

The  system  A. 10a  is  of  the  critical  type  if 

a)  at  rj  = 6 = 0 , B has  one  zero  eigenvalue 

b)  the  second  derivatives 


(0,0,) 


when  T]  = 6 = 0 


c) 


a / 0 . 


d2J 

oy1dy2 


vanish 


(A.  11) 


The  assumption  on  A. 10b  cai 
need  to  require  that  2 


dy- 


actually  be  weakened  slightly:  we  only 
vanish,  the  other  second  derivatives 


need  not  vanish.  However,  assumption  A. 10b  does  not  cause  any  loss 
of  generality  and  simplifies  the  analysis  considerably. 

In  terms  of  the  y coordinates,  for  small  rp 6 the  system 
(A. 10)  takes  the  form 

y2  = X(r),o)y2  + 0(y2,(ri  + 6)y) 

y1  = clti  + c26  + + c46)  (A.  12) 

+ y^Cgti  + c66)  + y1y2(c?n  + Cg'S’) 


A- 6 


(A. 12) 
cont ' d 


+ + ^io6)  + ayl  + ^yly2 


+ cyj^y 2 + 2fy  2 + 0(y4,(r|2  + 62)yL 


(ti  + 6)  y ) 


1 ^ 

Letting  y = — Y,  s=s/a  and  6^  = c2i-irl  + c2i^  we  ^ave 

y 2 = ^0y2  + °(y2'  <ti  + 6)y) 


y2  = Y1  + Y2yl  + Y3Y1  + Y4y.y; 


+ Y5y2  + yl  + byly2 


+ cyLy2  + dy2 


+ 0(y4,  (ti2  + 62)  y , 


(ti  + 6)y  ) 

We  now  introduce  new  variables  by 


(A. 13) 


Z2  = y2  + °(y  *(T)  + 6>y> 

Z1  = yl  + Slvyl  + fc  lr'yly2  + U1^2 

+ s26yl  + fc26yly2  + u26y2 

+ Wlyi  + W2yly2  + W3yly2 
+ w4y2  + 0 (y4 , (n2  + 62)y,riy3) 

2 

Without  loss  of  generality,  we  have  not  included  terms  0(y  ) in 
the  definition  of  Z^.  This  follows  from  assumption  A. 10b. 

A- 7 


(A. 14) 


From  A. 14,  we  have 


Z2  = y2  + 0 (yy  ' (rl  + 6)y  > 


zi  = yi  + 2si^Ly[  + t'liy1y2 


+ t,7iy[y2  + 2u,Tiy2y2  + 252yky^6 
+ t26(yly2  + yly2)  + 2u26y2y2 


+ 3wly^yi  + w2(ylyly2  + yly2> 


+ w3 (yly2  + 2yly2y2) 


+ 3w4y2y^ 


+ 0(y  3y,,Tl2y,,62y,) 


Using  (A. 13)  in  (A. 15)  yields  (with  X = X, 


Z2  = X0Z2  + 0(2  •r\z) 


z:  = 


Yi  + yx  y2 


+ y L ( y3  + 3wlYl) 


+ y1y2(Y4  + + t26^  + w2^i  + 2Xw3) 


+ y2<Y5  + 2ui  + 2u06\  + w-^X-t) 


3 1 


+ y. 


+ y^b  + w2\  + w2\  + v#2y2) 


A-8 


(A. 15) 


(A. 16) 


2 

+ y1y2(c  + 2w3^  + w3  y2) 


+ y2 (d  + 3w4\) 

+ 0 (y 4 , ( rj2  + 62)y,y3(n  + 6) 


We  choose 


W4  3X  w3  " 2\+Yo  W2  " X+y3 


Ulr'  + U26  = 


tlT1  + fc26  = 


w = ~ 
1 3 y i 


2X 

'Y4"W2Yl"2Xw3 


Noting  that  y3  = Z3  + 0(r|y3,n2y2)  and 


that  y^  = Z1y1  + 0(r]2Z2)  , egn  (A. 16) 
becomes 

Z2  = X0Z2  + 0 < z 2 » r]Z ) 

Z1  = Yi  + Y2Z!  + Z3  + 0 (Z4,  (r)2  + 62), 

(t|  + 6)  Z3)  , 


which  is  the  desired  normal  form. 
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